%30-nov-2007. This code reads Julia and LER & ESK data sets
% pick days with cont. common measurements
% 

%load c:\manoj\projects\plasma\HUA_PIU1.mat
load c:\manoj\projects\plasma\Julia_PLP.mat Julia_W fday w;
%load c:\manoj\projects\ace\ESKLER fday_esk_ler_ts LER_TS ESK_TS;
load c:\manoj\projects\ace\ESKLER fday_esk_ler_ts_filt FILT_LER FILT_ESK;
%esk_fday = fday_esk_ler_ts;
esk_fday = fday_esk_ler_ts_filt;
%FILT_ESK = (LER_TS-nanmean(LER_TS)+ESK_TS-nanmean(ESK_TS))./2;
FILT_ESK = (FILT_LER+FILT_ESK)/2;
%clear fday_esk_ler_ts LER_TS;
clear fday_esk_ler_ts_filt FILT_LER;
ESK_floor_fday = floor(esk_fday);
load c:\manoj\geomag\indices\aplist.mat;
mean_delta = (60/(24*3600));

%%

%Season Fliter
% L = season_filter(PN_fday,11,1,3,1);
% pcn_fday(~L) = [];
% PN_fday(~L) = [];
% PN(~L) = [];
% %season filter over---------


% Hd = chebi11_filter();
N_seg = 1;
N_data = 0;
mjd_date = datenum(2000,1,1);

for i = 1: 749,
    
    if sum(diff(Julia_W(i).k))+1== length(Julia_W(i).k) &...
        sum(isnan(w(i,:))) <= 60, 
    L = ESK_floor_fday == Julia_W(i).fday;
    if sum(L) > 0,
        
    esk_time = esk_fday(L);%+(4*5)/(60*24);%Try to advance 20 PCN time minutes
    esk_data = FILT_ESK(L);
    
%
    L = isnan(esk_data);
    if sum(L) < 10 & sum(L) ~= 0,
        pcn_data = interp1(esk_time(~L),esk_data(~L),esk_time);
    end;
%

    fprintf('Number of missiing points on fday %d = %d\n', Julia_W(i).fday, sum(L));
   %if abs(mean(diff(ace_time))-0.0028) <= 2.7778e-005, % Use this with Burke data
   %if abs(mean(diff(ace_time))-0.0035) <= 2.7778e-005, %Use this with ACE min averages
   if abs(mean(diff(esk_time)) - mean_delta) <= 2.7778e-005,% Use this with Thule 1 min PCN data
        JULI = w(i,Julia_W(i).k);
        ESKD = interp1(esk_time, esk_data, fday(i,Julia_W(i).k));
%         JULI1 = filter(Hd,JULI);
%         PCND1 = filter(Hd,PCND);
       
%         subplot(311);
        L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
        mean_ap = mean(ap(L));
        
   

        if sum(isnan(ESKD)) < 1,% & mean_ap >= 20,
            JULI_SEG(N_seg:N_seg+71) = JULI(1:72);
            ESK_SEG(N_seg:N_seg+71) = ESKD(1:72);
            TIME_SEG(N_seg:N_seg+71) = fday(i,Julia_W(i).k(1:72));
            N_seg = N_seg+72;
            N_data = N_data+1;
            data_length(N_data) = length(JULI);
            AP_SEG(N_data) = mean_ap;
 %           IMF_SEG(N_data) = mean_imf_bz;
               
         end;
   else,
       fprintf('Day %d has some missing time stamp\n', Julia_W(i).fday);
    end;
    end;
end;
end;

%[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(72),0,72,0.0033333);
[Cxx,F] = mscohere(JULI_SEG,ESK_SEG,hanning(72),0,72,0.0033333);
 
figure1 = figure;

axes('Parent',figure1,'XTick',[10 20 30 50 70 100 200 300 500 1000],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);

box('on');
hold('all');
xlabel('period in minutes');
ylabel('coherence');
hold on;
semilogx((1./(60*F(3:end))),Cxx(3:end),'r','LineWidth',2);

alpha = 0.99;
ci = 1- (1-alpha)^(1/(N_data-1));
aa = axis;
h=line([aa(1) aa(2)],[ci,ci],'LineStyle','-.','color','r');
axis([aa(1) aa(2) 0 1]);
title('(LER_H+ESK_H)/2  Vs Julia Vz');

    